Cluster Constraints on f(R) Gravity 
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Modified gravitational forces in models that seek to explain cosmic acceleration without dark 
energy typically predict deviations in the abundance of massive dark matter halos. We conduct the 
first, simulation calibrated, cluster abundance constraints on a modified gravity model, specifically 
the modified action f(R) model. The local cluster abundance, when combined with geometric and 
high redshift data from the cosmic microwave background, supernovae, Ho and baryon acoustic 
oscillations, improve previous constraints by nearly 4 orders of magnitude in the field amplitude. 
These limits correspond to a 2 order of magnitude improvement in the bounds on the range of the 
force modification from the several Gpc scale to the tens of Mpc scale. 



I. INTRODUCTION 

In /(i?) models, cosmic acceleration arises not from an 
exotic form of energy with negative pressure but from a 
modification of gravity. Here the Einstcin-Hilbert action 
is augmented by a general function f(R) of the Ricci 
or curvature scalar R [l], 0, [|[ . Such modifications not 
only can accelerate the background expansion but also 
generically lead to enhancements in gravitational forces 
on small scales. 

f{R) gravity is equivalent to a scalar-tensor theory, 
where fn = df jdR is the additional scalar degree of free- 
dom. This field has a mass and propagates on scales 
smaller than the associated Compton wavelength. Well 
within the Compton wavelength, the scalar mediates 
a 4/3 enhancement of gravitational forces, with corre- 
sponding strong effects on the growth of structure in 
the Universe. These enhancements are quantified by the 
mass of the field or equivalently by the value of the field 
in the background, fm- In order to pass stringent So- 
lar System constraints, viable f(R) models employ the 
chameleon effect which allows the field to become very 
massive in high-density environments However, in 
order for the chameleon effect to become active, the back- 
ground field should be smaller than the typical depth of 
cosmological potential wells, \fm\ < \^f\ ~ 1(U 6 — 1CU 5 

Independently of this Solar System constraint, how- 
ever, it is worth studying the constraints that can be 
placed on f(R) gravity directly from cosmological obser- 
vations. In the linear regime, these provide only weak 
constraints. Changes in the low multipole anisotropy of 
the cosmic microwave background (CMB) only place or- 
der unity bounds on the value of fua, while changes to 
the shape of the matter power spectrum, though larger, 
can be mimicked by galaxy bias Q . 

The effect of enhanced forces can be substantially more 
prominent in the non-linear regime. Cosmological sim- 
ulations have shown that for field values larger than 
\fna\ ~ 10~ 5 the abundance of rare massive halos are 



enhanced substantially 0]. Counts of galaxy clusters 
therefore provide the opportunity to improve cosmologi- 
cal constraints on f(R) models ultimately by 4-5 orders 
of magnitude. 

In this Paper, we quantify current cluster abundance 
constraints on f(R) models from a combined sample of 
low-redshift X-ray clusters. We begin in mi Al with a de- 
scription of the model and its impact on cluster predic- 
tions. In §1111 we describe the likelihood analysis of the 
local cluster abundance data. We combine these con- 
straints with geometric and high redshift constraints in 
§IVI to obtain upper limits on the modification to gravity. 
We discuss these results in §V] 



II. f(R) CLUSTER ABUNDANCE 

In this section, we describe the enhancement that f(R) 
models make to the cluster abundance. We begin in mi Al 
with a brief review of the f(R) model itself. We describe 
cosmological simulations with representative parameter 
choices in §11 Bl which are used to calibrate mass func- 
tion enhancements across a wider range of parameters in 
§HC1 



A. f(R) Model 

In the f(R) model, the Einstcin-Hilbert action is aug- 
mented with a general function of the scalar curvature 

R, 



R + f(R) 

16ttG 



(1) 



Here and throughout c = h = 1. Gravitational force 
enhancements are associated with an additional scalar 
degree of freedom fn = df / dR and have a range given by 
the Compton wavelength Ac = (Sdfii/dR) 1 ^ 2 . This ad- 
ditional attractive force leads to the enhancement in the 
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abundance of rare massive dark matter halos described 
below. 

For definiteness, we choose the functional form for 
f(R) given in [j| (with n = 1): 



f(R) = -2A 



R 



R + P 2 ' 



(2) 



with two free parameters, A, p 2 . Note that as R — > 
0, /(i?) — ► 0, and hence the model does not contain a 
cosmological constant. For R ^> p 2 , the function /(-R) 
can be approximated as 



f(R) = -2A-f R0 ^ 



(3) 



with /^o = ^2A^i 2 /^g replacing /i as the second param- 
eter of the model. Here we define -Ro = R{z = 0), so 
that fno = fn(Ro), where overbars denote the quantities 
of the background spacetimc. Note that if \fm\ <C 1 the 
curvature scales set by A = O(Ro) and p? differ widely 
and hence the R 3> p 2 approximation is valid today and 
for all times in the past. 

The background expansion history mimics ACDM with 
A as a true cosmological constant to order f R o- There- 
fore in the limit \fn,o\ "C 10~ 2 , the f(R) model and 
ACDM are essentially indistinguishable with geometric 
tests. Nonetheless geometric tests do constrain one of 
the two parameters (A) leaving the cluster abundance to 
constrain the other parameter (f R o) which controls the 
strength and range of the force modification. With the 
functional form of Eq. ([3]) , the comoving Compton wave- 
length becomes 



Ac 

1 + 2 

with a value today of 

Aco ~ 32 



Tro 



i? 3 



^4m P c 
10" 4 1 



(4) 



(5) 



The behavior of f(R) gravity is described by the mod- 
ified Einstein equations. Specifically, the trace of the 
linearized Einstein equations yields the f R field equation 



V 2 6fR = y [5R( f R ) - SnGSp m ] 



(0) 



where time derivatives have been neglected compared 
with spatial derivatives, coordinates are comoving, 5f R = 
Ir{R) - f R (R), 6R = R~R, S Pm = Pm - p m . Note that 
the local curvature R is given as a function of the lo- 
cal field value f R . The time-time component returns the 
modified Poisson equation 



V 2 * = AnGa 2 5 Pm - ^V 2 6f R . 



Here \1/ is the Newtonian potential or time-time metric 
perturbation 2^ = <5<?oo/.9oo m the longitudinal gauge. 



These two equations define a closed system for the New- 
tonian potential given the density field. The matter falls 
in the Newtonian potential as usual and so the modifica- 
tions to gravity are completely contained in the equation 
for 

Due to the conformal equivalence of f(R) and scalar- 
tensor theories and the conformal invariance of electro- 
magnetism, the geodesies of photons are unmodified by 
the presence of the scalar fg field save for conformal 
rescaling factors of 1 + f R Q. This means that given 
a fixed density field, e.g. a halo of mass M, the lens- 
ing potential will be identical to the one in GR in the 
\Ir\ ^ 1 limit that we work in. In other words, we will 
measure the "true" mass M through lensing. 

On the other hand, the mass inferred from dynamical 
measures, Md yn , such as velocities and virial tempera- 
tures is related to the dynamical potential which will be 
different in the presence of the f R field. In the low curva- 
ture limit where 5R -C 8nG5p m , Eq. ([7]) shows that the 
dynamical potential will be enhanced by 4/3. Hence the 
mismatch between Afd yn and M could be as large as 33%. 
Conversely, field fluctuations can saturate in deep grav- 
itational potentials as f R 0. Here 5R « 8irG6p m and 
force modifications disappear via the chameleon mecha- 
nism. Then, if the whole mass is contained in the satu- 
rated region, Md yn = M. We discuss the mass calibration 
of the cluster sample in Section IIIII 



B. Simulations 

We use a modified N-body simulation as described in 
@ and used in @, HU to quantify the impact of the force 
modification on the cluster abundance. Specifically we 
employ the system of equations defined by the modified 
Poisson equation ([7]) and the f R field equation ([6]) in the 
context of cosmological structure formation. The field 
equation for f R is solved on a fixed cubic grid, using 
a non-linear relaxation algorithm. The potential ^ is 
computed from the density and f R fields using the fast 
Fourier transform method. The dark matter particles are 
then moved according to the gradient of the computed 
potential, — V^. 

The simulations were performed for a range of back- 
ground field values \f R o\ = 10 -6 — 10~ 4 . We also simu- 
lated \f R o\ = which is equivalent to ACDM, using the 
same initial conditions. Note that the background expan- 
sion history for all runs is indistinguishable from ACDM 
to O(f R0 ). 

Cosmological parameters were fixed in the simulations 
to a flat ft A = 0.76, Q b = 0.04181, H = 73 km/s/Mpc 
model and initial power in curvature fluctuations A s = 
(4.82xl0~ 5 ) 2 at k = O^Mpc" 1 with a tilt ofn s = 0.958. 
All simulations possess 512 grid cells in each direction 
(7) and N p = 256 3 particles. Halos were identified in the 
simulations using a standard spherical overdensity halo 
finder 0]. In the next section, we describe our model 
for the f(R) effects on the halo mass function, which 
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allows us to extend predictions to a range of cosmological 
parameter values. 



C. Mass Function Enhancement 

In order to properly marginalize constraints over cos- 
mological parameters, we need a prediction of the mass 
function enhancement as a function of cosmological pa- 
rameters and the f(R) parameter /rq. Due to the large 
computing requirements for full f(R) simulations [9J, 
running simulations for a range of parameters is not 
feasible, and we use a model of the mass function en- 
hancement based on spherical collapse and the Sheth- 
Tormen (ST) prescription (see Q for details). Note that 
we use this prescription and the cosmological simulations 
that calibrate it to compute enhancements only. For the 
ACDM baseline predictions, we use mass function results 
from state-of-the-art numerical simulations (see mil Al) . 

The ST description for the comoving number density 
of halos per logarithmic interval in the virial mass M v is 
given by 



dn 



dlnMv 



M v jy J dlnAL 



where the peak threshold v = S c /a(M v ) and 



uf(u) = Asj^av 2 [l + {av 2 Y v ] exp[-au 2 /2] 



(8) 



(9) 



Here a(M) is the variance of the linear density field con- 
volved with a top hat of radius r that encloses M = 
47rr 3 p m /3 at the background density 



d 3 k 
(2nf 



\W(kr)\ 2 P L (k), 



(10) 



where Ph(k) is the linear power spectrum and W is the 
Fourier transform of the top hat window. The normal- 
ization constant A is chosen such that J dvf(y) = 1. The 
parameter values of p = 0.3, a = 0.75, and S c = 1.673 
for the spherical collapse threshold have previously been 
shown to match simulations of ACDM at the 10 — 20% 
level. The virial mass is defined as the mass enclosed at 
the virial radius r v , at which the average density is A v 
times the critical density p cv . For consistency with clus- 
ter analyses, overdensities will refer to critical density 
throughout the paper; the corresponding overdensities 
in terms of the background matter density are given by 

A Pm = A Pc jn m . 

Ref. 0] derived a model for the mass function enhance- 
ment measured in the f(R) N-body simulations. The 
mass function calculation is based on the Sheth-Tormen 
prescription using the linear power spectrum for the f(R) 
model, and two limiting cases for the spherical collapse 
parameters. In one case, we simply assume that the 
spherical perturbation considered is always larger than 
the (local) Compton wavelength of the fa field, so that 



Full simulation, |f B0 |=10- 4 

Spherical collapse 
Mod. force spherical collapse 




FIG. 1: Mass function enhancement at z = with respect 
to ACDM as a function of M = M500, measured in f(R) 
simulations with |/i?o| = 10 -4 . Also shown is the range of 
spherical collapse predictions from Q. For the constraints, we 
conservatively use the lower limit of the shaded band (dashed 
line) . 



gravity is GR throughout, and the spherical collapse pa- 
rameters are unchanged (5 C = 1.673 and A v = 94 for 
£l m = 0.24). In the second case, we assume that the 
perturbation is always smaller than the local Compton 
wavelength, so that forces are enhanced by 4/3. The 
corresponding modified spherical collapse parameters are 
S c = 1.692 and A v = 74 for Q m = 0.24. Fig. Q] shows 
the range of predicted mass function enhancement for 
these two limiting cases, and the results of the f{R) sim- 
ulations, for |/flo| = 10~ 4 . The mass definition used, 
M = Ma with A = 500 relative to the critical density, is 
the same as used in the X-ray cluster measurements. In 
order to obtain conservative upper limits, we choose the 
modified force prediction for the mass function, which 
corresponds to the lower bound of the shaded band in 

Fig.m 

For a given set of cosmological parameters (A s , fl m , h, 
/bo)) we nrs t calculate the Sheth-Tormen mass functions 
dn/dlnM v for ACDM and f(R) using the respective lin- 
ear power spectra. We then rescale each mass function 
from the respective virial mass to the common mass def- 
inition M = -M500, using the procedure outlined in [ill ]. 
We need to assume a halo profile for this mass rescaling, 
which we take to be of the NFW form with the concentra- 
tion relation given by [l§]. As shown in Q, the profiles 
of dark matter halos in f(R) within r v are sufficiently 
similar to those measured in GR simulations that f(R) 
effects can be neglected in the mass rescaling. Finally, A v 
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for ACDM is obtained from the fitting formula of [l3[ : 

AG R (fi m ) = 178f7° m 45 . (11) 

For the f(R) enhanced forces, we assume the same scal- 
ing, fixing the ratio Al iR) /A^ R = 74/94. We neglect the 
small f2 m dependence of the collapse thresholds within 
the range of interest, and keep <5 C (ACDM) = 1.673 and 
6 c {f(R)) = 1.692 fixed. 

III. CLUSTER LIKELIHOOD 

In this section we describe the cluster likelihood as a 
function of cosmological and f{R) parameters. Since we 
assume a spatially flat cosmology and the expansion his- 
tory of f(R) models are indistinguishable from ACDM, 
the main cosmological parameters that we have to con- 
sider are J7 m , h, and the primordial normalization A s at 
k = 0.02 Mpc -1 . Other parameters such as the power 
spectrum tilt do not affect the constraints appreciably 
jl4| . Since the CMB determines f2 m /i 2 to good preci- 
sion, we are mainly dealing with f2 m , A s , and the f(R) 
parameter fnQ for the cluster abundance. 

In 9111 A[ we review the likelihood approach taken in 
Ref. [l4[. In i jlllBl we describe the rescaling of these 
results for f(R) models. 

A. ACDM 

The cluster sample used in this work is the low- z sub- 
sample of 49 clusters described in 0, [l5|. This is an 
X-ray flux-limited sample of clusters originally detected 
in the ROSAT All-Sky Survey at high Galactic latitudes 
and at z > 0.025. All of the objects were later observed 
with Chandra, providing low-scatter proxies for the total 
mass which can be constructed from the mean X-ray tem- 
perature and gas mass (see below). The effective rcdshift 
depth of this sample is z < 0.15. 

Cluster total masses are estimated using the Yx pa- 
rameter defined as Yx = Tx x M gas , where Tx is the 
average temperature measured from the integral X-ray 
spectrum, and M gas is the estimated gas mass derived 
from the analysis of the X-ray surface brightness profile. 
Yx is a direct X-ray observable, even though it is hard to 
backtrack it to raw observablcs such as the total X-ray 
luminosity. For a detailed description of the data analysis 
procedures, see fl5| . 

Based on state-of-the-art cosmological numerical sim- 
ulations, Yx is expected to be tightly correlated with 
the total cluster mass, M cx Y£ /S H(z)~ 2 / s , with < 
10% scatter Numerical experiments show that the 
Yx — M relation is remarkably insensitive to the cluster 
dynamical state. The power law slope and evolution fac- 
tor are also insensitive to the details of star formation 
history and non-gravitational heating of the intracluster 
medium although these processes do change the overall 



normalization of the relation (e.g, [l6|, [13] )■ The normal- 
ization of the Yx — M relation was determined observa- 
tionally [l5| using two independent techniques, (1) by X- 
ray hydrostatic method using a subsample of dynamically 
relaxed clusters, and (2) by weak lensing mass measure- 
ments. The two methods are in excellent agreement, and 
this was used to place an upper limit on systematic un- 
certainties in the cluster mass scale, AM/M < 9%. For 
our purposes, this agreement also means that the normal- 
ization of the Yx — M relation is tied to the weak lensing 
measurements, which should provide the true mass in the 
f(R) theories we consider (§ 111 A[) . 

The cluster component of the likelihood function is 
computed assuming a purely Poissonian nature of sta- 
tistical fluctuations 1 : 

lni = lnp(Mf * , Zi) + ln M i St 

, , 1 (12) 
- dz dM cst p(M cst ,z), 

where M? st and is the estimated mass and rcdshift of 
cluster i, p(Mf st : Zi) is the model probability density to 
observe a cluster with mass M° st at redshift Zi , and sum- 
mation is over the clusters in the sample and integration 
is over pre-selected z m - m — z max and M m - ln — M max in- 
tervals. Note that the ln Mf st term appears because the 
mass estimates and hence the mass binning required to 
convert probability densities into probabilities depends 
on cosmology. 

The model probability density, p, is a product of the 
Tinker ct al. mass function model [Tc| for a given set 
of ACDM parameters, cosmological dV/dz function, and 
the survey selection probability as a function of mass 
and redshift. The product of these components is con- 
volved with the intrinsic and measurement scatters in the 
Yx — M relation. The computation of all these terms is 
discussed in detail in [l5[ , and the ACDM parameter con- 
straints derived from this cluster dataset are presented in 
[l4| . We now turn to the computation of the cluster like- 
lihood in the f(R) models. 



B. f(R) Scaling 

A full likelihood analysis of the f(R) cluster constraints 
would entail a recalculation of the X-ray cluster likeli- 
hood grid [HI . Since the modification to the shape of the 
cluster mass function in f(R) is not dramatic across the 
dynamic range of masses probed by the cluster mass func- 
tion data (e.g., Fig. 16, 17 in [TB]), we opt for a simpler 
approach: for each point in parameter space Q m , A s , fno, 
we calculate the f(R) mass function enhancement at a 



1 We ignore a contribution from cosmic variance; the validity of 
this assumption is justified in Il5ll . 
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FIG. 2: Mass function enhancement for \fno \ = 10~ 4 from the 
spherical collapse model (black, solid) as in Fig. [TJ and the 
corresponding enhancement when increasing the linear power 
spectrum normalization in ACDM. The vertical line indicates 
the pivot mass M E g used to calculate the likelihood. The blue 
dashed line shows the enhancement for a rescaled normaliza- 
tion (cr§ ff — as x 1.066) that matches the f(R) enhancement 
at M eff . 



FIG. 3: Cluster likelihood as function of primordial normal- 
ization (quantified by the linear power spectrum normaliza- 
tion cr^ CDM a ACDM model would give), for fixed values of 
On = 0.258, h = 0.716, and \ fm \ = 10 -4 in case of the f(R) 
prediction. The red short-dashed line shows the likelihood 
calculated using the full f(R) mass function enhancement, 
while the blue long-dashed line shows the ACDM likelihood 
obtained with the rescaled normalization, erf . 



pivot mass, M cS = 3.677 x 1CT 14 M Q /h. Here we take 
M to be the true or lensing mass, which is the most 
conservative assumption. Equating the dynamical mass 
Mdy n ~ 4Af/3 to the Yx based mass can only increase 
the abundance enhancement in the f{R) models. We cor- 
respondingly also ignore the additional tension in f(R) 
implied by the observed agreement between the lensing 
and Yx masses. 

Then, for this set of parameters, the enhancement is 
converted to an effective, not actual, linear power spec- 
trum normalization, er| , assuming a ACDM model with 
the same Q m and the mass function prescription of [l8[. 
This approximation assumes that, in the mass range 
probed by the X-ray clusters, the mass-dependence of 
the mass function enhancement due to f(R) has the same 
shape as that due to an increased power spectrum nor- 
malization. This is only approximately true (see Fig. [5]), 
since the growth is scale-dependent in f(R). 

In order to benchmark the accuracy of this simplified 
approach, we recalculated the cluster likelihood for fixed 
values of Q m ,h and a range of A s using the full f{R) 
mass function, for |/ro| = 10~ 4 . We then compared this 
likelihood with the ACDM cluster likelihood calculated 
for the rescaled normalization, cr§ . The resulting like- 
lihoods are shown as function of the primordial normal- 
ization in Fig. [3J First, clusters clearly prefer a lower 
primordial normalization in the f(R) model, due to the 



enhanced growth. The approximate likelihood calculated 
using the rescaled normalization agrees quite well with 
the full likelihood. Note that this approach is in any case 
conservative, as the constraints are weakened (though not 
significantly) by neglecting the additional information. 

While all f(R) mass function enhancements were cal- 
culated at z = 0, we verified that the evolution of en- 
hancements in the rcdshift range probed by the cluster 
sample, z = — 0.15, is negligible. 

Fig.[3]shows that if S7 m w 0.26 and the primordial nor- 
malization of the simulations were verified to very high 
precision by external constraints, and the mass calibra- 
tion of the cluster sample carried no systematic error, 
the cluster abundance would be able to rule out the f(R) 
model with \fm\ = 10~ 4 at around 95% confidence. We 
next address to what extent these expectations apply to 
the joint cosmological and cluster data. 



IV. COMBINED CONSTRAINTS 

The excess cluster abundance predicted in f(R) mod- 
els can be converted into limits on the field amplitude 
fno once data external to the clusters has fixed the back- 
ground expansion history and primordial normalization 
of density fluctuations. In mV Al we describe the external 



6 



data that we use for these purposes and present combined 
results in gVBl 



A. External Data Sets 

CMB: Following Ref. [3], we employ a simplified ap- 
proach to incorporating CMB constraints from WMAP5 
into the cluster analysis. We take three CMB parame- 
ters — angular scale of the first acoustic peak, £a', the 
so called shift parameter, R; and the recombination red- 
shift, z*. The likelihood for the geometric side of the 
WMAP-5 data is computed using the covariance matrix 
for I a, R, and z* provided in [lj|. In the ACDM ex- 
pansion history these quantities depend on O m , h and 
fl b h 2 . 

Next we add the information contained on the initial 
amplitude of fluctuations. The WMAP team provides 
the amplitude of the curvature perturbations at the k = 
0.02 Mpc -1 scale, 



A s = (2.21 ± 0.09) x 10" 



(13) 



To implement this constraint in terms of er^ CDM and the 
chosen cosmological parameters we use the fitting for- 
mula 0: 



a\' 2 



T ACDM 



1.79 x 10 4 V 0-024 
x (7.808 /i) (1 "™ )/2 



1/3 



0.14 

0.693 



0.72 



-0.563 



0.76 



(14) 



(we adjusted numerical coefficients to take into account 
that the original fit uses the CMB amplitude at k = 
0.05 Mpc -1 while the WMAP-5 results are reported for 
k = 0.02 Mpc -1 ). In this equation, Go is the growth 
suppression relative to S oc (1 + z)^ 1 due to A evaluated 
today. We then include a x 2 contribution given by 



XCMBnorm = (A. X 10 9 - 2.21) 2 /0.09 2 



(15) 



The XcMBnorm component is computed assuming a fixed 
n = 0.95; the results are completely insensitive to varia- 
tions of n within the WMAP measurement uncertainties 
and even to setting n = 1 . The sum of the geometric and 
growth component of the CMB % 2 is marginalized over 
Qbh 2 and h. The end result is a Xcmb f° r the CMB that 
is a function of f2 m and cr^ CDM . 

SN: We use the distance moduli estimated for the 
Type la supernovae from the HST sample [2l[ , SNLS sur- 
vey HU, and ESSENCE survey [H], combined with the 
nearby supernova sample as compiled in Ref. [2~i| . Calcu- 
lation of the SN la component of the likelihood function 
for the given cosmological model is standard and can be 
found in any of the above references. The end result is a 
Xg N that depends on Q. m only. 

Hq: We use the recent determination of Hq [25| . 
H a = 74.2±4.8km/s/Mpc, in conjunction with the CMB 



Clusters+CMB 




FIG. 4: Top panel: Likelihood contours from clusters+CMB 
in the fm — fl m plane, marginalized over the primordial 
normalization. Shown are 68.3%, 95.4% and 99.7% likeli- 
hood contours. Bottom panel: Likelihood contours in the 
fm) — f2 m plane marginalized over the primordial normaliza- 
tion, for clusters+CMB only and when including other geo- 
metric measurements. 



constraint of Q, m h 2 = 0.133±0.006 \M as a measurement 
of f2 m . Marginalizing over the uncertainty in Q, m h 2 re- 
sults in the following Gaussian likelihood: 



Xh 



n m - 0.242 
0.034 



(16) 



BAO: In a similar way, we use the distance scale given 
by the baryon acoustic oscillation measurements (BAO) 
of [27[. We use their Eq. (16), which after marginalizing 
over Vl m h 2 yields: 



Xbao 



0.285 



0.019 



(17) 



The BAO constraint is in fact the most precise one and 
hence dominates our O m likelihood. 

Finally, we combine all the contributions of external 
data sets 

xLt = Xcmb + Xh + x!n + xIao . ( 18 ) 
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CMB+clusters 
CMB+clusters + SN + H n 



0.001 



0.002 



0.003 



FIG. 5: Likelihood relative to f m = (ACDM) as a func- 
tion of fm for CMB+clusters and including other geometric 
measures. We have marginalized over fi m and the primordial 
normalization. The horizontal lines show the 2a and 3a con- 
fidence levels. Using all measures combined, |/ho|/10 -4 < 1.3 
at 95% confidence level. 



and add Xext to the cluster likelihood contribution of 



-21nL. 



B. Results 

In Fig. 2] (top) we show the results of combining the 
cluster abundance and CMB constraints. The assump- 
tion of spatial flatness in combination with the CMB data 
alone constrains fi m and limits the extent of the /#o — fi m 
degeneracy. Note that the bounds on fm tighten as O m 
increases. With only these two data sets the statistical 
upper limit after marginalizing over O ro is |/ro|/10 -4 < 8 
at 95% statistical confidence level (CL). 

Data on SN distances, Hq and BAO distances tighten 
the bounds on £l m reducing the degeneracy with |/i?o|- 
Fig. [4] (bottom) shows that the BAO data in particular 
make a strong impact on constraints since they favor high 
fl m . In Fig. [5]we show the Sy 2 statistic after marginal- 
ization of Q, m . With all of the data, |/ro|/10~ 4 < 1.3 at 
95% statistical CL. Table fl] summarizes the upper limits 
on fno for the different data sets and for different confi- 
dence levels. 

The main caveat to these statistical constraints is the 
possibility of systematic shifts in the mass calibration 
of the cluster sample. In Fig. [5] we show the impact 
of ±9% shifts in the cluster mass scale on the clus- 
ter+CMB constraints. Note that this effect mainly shifts 



10- 



10" ; 




Clusters+CMB 
+9% mass shift 
-9% mass shift 



0.28 



0.3 



FIG. 6: Likelihood contours from clusters+CMB in the 
fno — H m plane, marginalized over the primordial normaliza- 
tion as shown in Fig. 3] The solid contours show the results 
for the standard cluster mass scale, while the black solid (red 
dashed) lines show the likelihood in case cluster masses are 
overestimated (underestimated) by 9%. 



the contours by Af2 m ~ ±0.015. If we assume that clus- 
ter masses are underestimated ("—9% mass shift", i.e. 
Mx,obs = 0.91 M), the abundance at a fixed mass is in 
fact higher, and hence allows higher /ro values. Con- 
versely, in the case that cluster masses are overestimated 
("+9% mass shift", i.e. M x ,obs = 1.09 A/), the true 
abundance at a fixed mass is smaller, hence tightening 
fno constraints. 

We show the impact of a ±9% mass calibration error 
on the final joint results in Fig. FjJ and in Tab. [U The net 
result is that the 95% statistical CL carries systematic 
errors of +1.7/ — 0.6 x 10 -4 , which we shall write as 
\f RO \/lO- 4 < 1.3±l*. 

Note that given M& yn m 4M/3 in f(R) if there is no 
screening due to the chameleon mechanism, the X-ray 
measurements possibly overestimate cluster masses by up 
to 33%. Hence, we expect that our use of lensing masses 
in calibrating the enhancement makes our constraint con- 
servative even given the possibility of a 9% underestimate 
in the X-ray mass measurement. 

Our model of the mass function enhancements (i fll C|) 
also represents a lower bound which always underpre- 
dicts the enhancement measured in N-body simulations 
for 10~ 6 < | /ho | < 0]. We have also determined 

upper limits on \fno\ using the less conservative limit- 
ing case presented in Q, which corresponds to using al- 
ternate collapse parameters that correspond to the GR 
values of S c and A„. This case is shown as the upper 
boundary of the shaded band in Fig. [2J While the pre- 
diction is still below the simulations, and a better fit, 
for | /ho | ^ 10~ 4 , it overpredicts the enhancements for 
smaller field values 0- The last line in Tab. Q] shows the 
resulting limits, which are tighter by a factor of 3-4. We 
cannot guarantee that this model does not ovcrpredict 
the enhancement in some region of the parameter space 
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0.001 0.002 0.003 

|f R ol 

FIG. 7: Likelihood as a function of fjio as in Fig. [S] for 
CMB+clusters and all geometric measures. The shaded 
band shows the weakening/strengthening of the constraints 
when varying the absolute mass scale by ±9%, correspond- 
ing to the estimated systematic uncertainty. In the weak- 
ened case (masses underestimated), the constraint degrades 
to I/mI/KT 4 < 3.0 at 95% confidence level. 

involved, but the corresponding tightening of the con- 
straints again indicates that our quoted limit should be 
viewed as conservative. 



V. DISCUSSION 

We have provided the first, simulation calibrated, clus- 
ter abundance constraints on a modified gravity model, 
specifically f(R) gravity. Enhanced forces below the 
Compton wavelength of the scalar field lead to corre- 
sponding enhancements in the cluster abundance, mak- 
ing the latter a sensitive probe of gravity on cosmological 
scales. 



TABLE I: Upper limits on fno in units of 10 



Confidence level (CL) 


68.3% 


95.4% 


99.7% 


Clusters+CMB 


1.0 


7.9 


> 31 


Clusters+CMB+SN+tfo 


1.0 


5.4 


> 31 


Clusters+CMB+SN+i/o+BAO 


0.3 


1.3 


6.3 


with +9% mass shift 


0.2 


0.7 


3.1 


with —9% mass shift 


0.7 


3.0 


14.7 


with alternate collapse parameters 


0.09 


0.4 


1.8 



Combined with constraints on the primordial ampli- 
tude of fluctuations from the CMB, and geometric con- 
straints from CMB, supernovae, Hq and BAO, the clus- 
ter abundance provides the most stringent constraint on 
these enhancements to date. In terms of the field ampli- 
tude today, the constraint is |/ro|/10~ 4 < l.S^g, where 
the range reflects a =f9% mass calibration error, an im- 
provement over previous constraints by 4 orders of magni- 
tude. This corresponds to an upper limit on the range of 
the gravitational force modifications of Ac < 38^® Mpc 
in this f(R) model. 

Our constraint should be viewed as conservative even 
given the 9% mass calibration error. We have ignored 
an overestimate of dynamically based A-ray masses over 
true or lensing masses, which could be up to a 33% 
shift that would further enhance the abundance and 
strengthen the constraint. In addition, we have not con- 
sidered the possibility of constraining f(R) force modi- 
fications from the comparison of cluster lensing and dy- 
namical masses. 

Furthermore, our model of the mass function enhance- 
ments represents a conservative lower bound which al- 
ways underpredicts the enhancement measured in N- 
body simulations. We found that the constraints tighten 
significantly if we use a less conservative model but a ro- 
bust implementation would require more accurate mass 
function calibration across the parameter space. 

On the observational side, current constraints are lim- 
ited mainly by systematics in the mass calibration and 
secondarily by the small number of local clusters. Relax- 
ing the assumption of a flat universe is not expected to 
degrade the upper limits appreciably. This is because our 
constraints only depend strongly on the allowed range in 
fl m and marginalizing over curvature changes this range 
negligibly once BAO are combined with SN and/or the 
CMB (23|. 

In the future, the abundance of massive clusters could 
ultimately provide another order of magnitude improve- 
ment in the field amplitude to |/ro| ~ 1CP 5 , rivaling 
solar system tests of gravity but in a very different, low 
curvature regime Q. Below field values of ~ 10 -5 , the 
chameleon mechanism suppresses the enhancement at 
the high mass end 0]. In this regime, further improve- 
ments are potentially available if the abundance of galaxy 
groups can provide constraints on the halo abundance at 
intermediate masses. 

While we have considered a specific functional form 
of f(R) here @, different functional forms have been 
proposed in the literature, see e.g. [H, HH, H3|. These 
models differ primarily in the evolution with rcdshift of 
the Compton wavelength of the fn field. Hence, we ex- 
pect our results to be generic once the field amplitude 
and range are rescaled to some effective redshift which 
matches the impact on the linear growth today on a scale 
relevant for clusters. For example in models where the 
curvature dependence is steeper, the field amplitude to- 
day is allowed to be larger since its current value has 
little impact on the growth of structure. In these models 
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solar system tests become even more powerful relative to 
cosmological tests. 

More generally, the abundance of galaxy clusters 
promises to be a good probe of other modified gravity 
scenarios as well, such as DGP and other braneworld 
models once their mass functions are calibrated by cos- 
mological simulations [3l|, . 
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